Astron. Nachr./AN xxx (xxxx) x, xxx-xxx 



Differential rotation and meridional flow in the solar convection 
zone and beneath 



L.L. KiTCHATiNOV^'2 and G. Rudiger^ 



in 
o 
o 

(N 



^ Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482, Potsdam, Germany 
^ Institute for Solar-Terrestrial Physics, PO Box 4026, Irkutsk 664033, Russia 



Received 30 April 2005; accepted 30 May 2005; published online 1 July 2005 



> 

cn 
\o 

O 

in 
o 

o 



Abstract. The influence of the basic rotation on anisotropic and inhomogeneous turbulence is discussed in the context of 
differential rotation theory. An improved representation for the original turbulence leads to a yl-effect which complies with 
the results of 3D numerical simulations. The resulting rotation law and meridional flow agree well with both the surface 
observations (dO/dr < and meridional flow towards the poles) and with the findings of helioseismology. The computed 
equatorward flow at the bottom of convection zone has an amplitude of about 10 m/s and may be significant for the solar 
dynamo. The depth of the meridional flow penetration into the radiative zone is proportional to fco?c! where I'coro is the 
viscosity beneath the convection zone. The penetration is very small if the tachocline is laminar. 
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1. Motivation 

The internal solar rotation is well-known from helioseismo- 
logy (Wilson et al. 119971 Kosovichev et al. 1997; Schou et 
al. .1998 1. The decrease of angular velocity with latitude ob- 
served on the solar surface survives throughout the convec- 
tion zone. The differential rotation decreases sharply with 
depth in the thin 'tachocline' beneath the convection zone. 
The helioseismology inversions also show a considerable in- 
crease of rotation rate with depth just beneath the solar sur- 
face (Fig.Q. 

The present paper demonstrates how mean-field models 
work to explain the overall pattern of the internal solar rota- 
tion and to derive the related meridional flow. To this end 

- our former theory of angular momentum transport by ro- 
tating turbulence (Kitchatinov & Rudiger[l993j hereafter 
KR93) is extended to the case of slow rotation 

- the global stellar circulation model is modified to in- 
clude the meridional flow penetration below the convec- 
tion zone. 

The currently debated penetration (Nandy & Choudhury 
120021 Oilman & Miesch I2004> may be important for 
advection-dominated models of the solar dynamo (Choud- 
huri, Schussler & Dikpati 1995 Dikpati & Oilman 2001 ; Bo- 
nanno et al. 2002 1. The dependence of the penetration depth 
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on basic parameters is computed with our global stellar cir- 
culation model. The new simulations include also a computa- 
tion of the solar tachocline with one of the currently discussed 



tachocline models (Riidiger & Kitchatinov l 1 9971 MacOregor 
& Charbonneau 1999/ 

The angular momentum fluxes inside the convection zone 
basically depend on the Coriolis number, 

n* = 2reorr/?, (1) 

where J? is the angular velocity and Tcon is the convective 
turnover time. The numerical models (Kitchatinov & Riidiger 
[T^ [1999; Kuker & Stix 2001) based on the A-effect theory 
of KR93 had practically no free parameters but reproduced 
the main features of the solar rotation law and predicted cor- 
rectly the trend of the surface differential rotation with stellar 



mass and rotation rate (see Hall 1991; Barnes et al. l2005t . 

The characteristic value of the Coriolis number for the 
Sun exceeds unity, ~ 6, and it is even larger for rapidly 
rotating stars. The Coriolis number varies, however, with 
depth inside the Sun (cf. Fig. |5J and becomes rather small 
near the surface. The former models did not concern the near- 
surface region (e.g. Kitchatinov & Riidiger' 1995 1. 

The present paper extends the former results by using a 
more general representation for the spectral tensor of turbu- 
lent convection. The resulting expressions for the yl-effect in- 
clude a free parameter which is unimportant for fl* ^ 1. Its 
value can be restricted by numerical experiments. The new 
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Fig. 1. Internal solar rotation from helioseismology. Courtesy 
NSF's National Solar Observatory. The numbers give the lat- 
itude rather than the colatitude. 
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Fig. 2. Profile of the Coriolis number Q in the solar model 
ofStix&Skalev (.1990j . 



model shows agreement with both the observed rotation law 
and the meridional flow near the solar surface. It predicts the 
meridional flow towards equator with an amplitude of about 
10 m/s at the bottom of the convection zone. 



2. Angular momentum transport 
2.1. Basic relations 

The ability of turbulence to transport angular momentum is 
reflected in the structure of the correlation tensor, 



Q^ 



{u[{x, t)u'Ax, t)), 



(2) 



of the fluctuating velocity u'. The radial and the latitudi- 
nal angular momentum fluxes are proportional to the off- 
diagonal components, Qr^, Qotp, of the tensor in spherical 
coordinates. The fluxes are finite even for the case of rigid 
rotation; they are conventionally parameterized as 



The eddy viscosity, v^, is introduced here for dimensional 
reasons, V and H are the normalized vertical and horizontal 
fluxes. 

Negative V also for slow rotation is required to explain 
the negative radial gradient of the angular velocity seen in 
the subsurface layer in Fig. ^ The equatorial acceleration, 
on the other hand, demands positive H (Riidiger 1989). For 
a nonuniform rotation the correlation tensor (|2j includes the 
diffusive part, , which is defined by the eddy viscosity isn- 
sor,7Vyfei, 

Q^J=QtJ+Q''^r (4) 

with 

Qr, = -M...£^, (5) 

where u is the mean velocity. 

2.2. 3D simulations 

It is not easy to separate the yl-effect from its viscous coun- 
terpart in the results of 3D numerical simulations. Kapyla et 
al. ( i2004t suggest that the diffusive momentum fluxes may be 
small in the box simulations. They are also reporting H > Q. 

Already PuUckinen et al. (1993) found that the function 
V is negative and H is positive. The results of Rieutord et 
al. (1994) confirmed the positive H but the vertical flux V 
was small. Chan (2001) found very clear results also with 
positive H but with negative V. The equatorial value of the 
latter proved to be rather small. 

Simulations of compressible thermal convection under 
the influence of rotation were also made with the finite- 
difference, fractional-step code NIRVANA (version II) in a 
small rectangular box defined on a Cartesian grid (Ziegler 
1998, 1999; Rudiger et al. l2005b> . The Coriolis number Q 
varied from about unity at the top of the convection box to 
about four at its base. The results agree with those of Chan 
(2001 ) and Kapyla et al. (20041 in providing negative V and 
positive H. The horizontal flux is strongly concentrated to- 
wards the equator 

The box simulations can be summarized as follows. The 
radial flux V of the angular momentum is always negative. 
Its equatorial value vanishes for SI* > 1 in agreement with 
KR93 for fast rotation. The horizontal flux H is positive in 
all simulations. It vanishes at the poles as it should be from 
symmetry reasons. 

2.3. Quasilinear theory 

The quasilinear theory of the angular momentum transport 
by rotating turbulence in density-stratified fluids results in the 
expressions 



(6) 



for the normalized fluxes'. The coefficients V^'^*'' and H^^\ 

2 



^(1) 



{Jo{n*) + aio{n*))., 



{Ji{n*)+ah{^*)). 



(7) 



Qed) — i^T^ H cos ( 



(3) 



see the Appendix for more details 



3 THE MODEL 
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Fig. 3. The functions of the Eq. 0. Only Jq and Iq remain 
finite at small Coriolis numbers and they have opposite signs. 
For rapid rotation, i.e. large Ji dominates all other func- 
tions. 
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Fig. 4. Normalized fluxes of angular momentum (|6j in de- 
pendence on latitude for a = 2 and £7* from 0.3 to 10. The 
lines are marked by the correspondent values of the Coriolis 
number 

depend on the Coriolis number via the functions J and / 
given by ( IA18> and (A19). The functions Jq and J\ are 
the same as in KR93. All these 'v4-functions' are plotted on 
Fig-E It can be seen from the plots and Eq. Q that a positive 
i?*^^^ requires a slightly positive anisotropy parameter a de- 
fined by ( IA12> and ( IA15> . On the other hand, a negative V for 
slow rotation can be reproduced only with a > 1. The van- 
ishing V at the equator for rapid rotation is due to the rapid 
decrease of the functions Jo and /o with il* . 

Figurel^shows the normalized fluxes, V and iJ, as func- 
tions of latitude for a = 2,a value that was also used to model 
differential rotation and meridional flow as discussed below. 
The resulting profiles are in agreement with all the main find- 
ings of the box simulations, i.e. 

- positive H, 

- increasingly negative V with increasing latitude, 

- negative equatorial V approaching zero with increasing 
rotation rate. 



The anisotropy relation 



follows from (A16) which for a = yields positive (nega- 
tive) values in the upper (lower) part of the solar convection 
zone. With a ~ 2 the amplitude of this anisotropy is strongly 
reduced in the upper part of the convection zone. In its lower 
part, the effect is completely quenched by the rapid rotation 
and the turbulence there remains vertically-dominated. 

3. The model 

The global flow in the solar convection zone was computed 
with our former model (Kitchatinov & Riidiger 1999 1 which 
now applies the new expressions Q for the A-effect with 
the functions ( lAlSt and (A19). The model solves the steady 
mean-field momentum equation, 

{u-\/)u+-\/ ■{pQ) + -\/P = g, (9) 
P P 

for axisymmetric mean flow u simultaneously with the en- 
tropy equation 

I T-irad I 



pTu ■ \/S + div F 



0. 



(10) 



The entropy equation is needed to define the baroclinic part of 
the meridional flow as well as the eddy transport coefficients 
for the nonrotating case, e.g.. 



15Cr, 



dr 



(11) 



The model accounts for the anisotropy and quenching of the 
eddy heat transport coefficients due to global rotation (see 
Riidiger et al. 2005a, for details). 

The usual stress-free and zero-penetration boundary con- 
ditions are applied. The lower boundary is located at the base 
of the convection zone (Model 1). In this model the boundary 
conditions exclude any penetration of meridional flow in the 
stably stratified layer beneath the convection zone. The mag- 



2^2 , 
corr ^ 



(8) 



netic tachocline model by Riidiger & Kitchatinov J1997> is 
also used in Model 1 so that the rotation law is computed for 
the entire volume of the Sun. The model realizes the magnetic 
tachocUne as a Hartmann layer due to a weak internal mag- 
netic field of the solar radiative core. The tachocline mod- 
eling has no influence on the rotation law computed for the 
convection zone proper. It only uses the results of these com- 
putations as the top-boundary condition for the tachocline. 

Penetration of the meridional flow into the stably strati- 
fied layer beneath the convection zone has been probed in 
former models but it was always weak. We now rediscuss 
this issue because the penetration became important in re- 
lation with new ideas of solar dynamo theory. The bottom 
boundary is thus placed beneath the convection zone inside 
the stably stratified region (Ekman layer model. Model 2). 
For numerical reasons, an extra but small isotropic and uni- 
form viscosity and thermal conductivity were prescribed for 
the entire computation region. The code is flexible enough to 
apply a nonuniform grid with about 300 grid points beneath 
the convection zone to resolve even a shallow penetration. 
The dependence of the penetration depth, -Dpcn, on the pre- 
scribed effective viscosity (i^coic) is computed with Model 2. 
The computation of a tachocline makes little sense in combi- 
nation with this model and it was not attempted. 



4 RESULTS AND DISCUSSION 



4. 1 Global How structure 



With the thermal boundary conditions both spherically 
symmetric heating from below at the bottom of the compu- 
tational domain and black-body radiation of the photosphere 
are introduced into the model. The top boundary of the model 
is placed shortly below the photosphere to exclude the sur- 
face layer with very sharp stratification for which case our 
representation for the A-effect is certainly incomplete. Nev- 
ertheless, in order to include the supergranulation layer into 
the simulations, the external boundary was shifted from its 
usual position at — 0.95 to — 0.98. 

Yet, the surface layers are now unstable to large-scale 
thermal convection with the usual value ^corr/^p — 1-7 of 
the mixing length auLT (Tuominen et al. ll994> . We had to 
increase the parameter to qmlt = 2.1 to avoid this instabi- 
lity. With such a large value the mixing length is larger than 
the density scale height, which is problematic for the con- 
struction procedure of the spectral tensor used in the Ap- 
pendix. 

The solar model of Stix and Skaley ( I1990> is used. The 
temperature and density at the external boundary are po = 
1.37 • IQ-^ g cm-3 and = 1.03 • 10^ K, taken as the input 
parameters of the simulations. 

4. Results and discussion 
4.1. Global flow structure 

After a confrontation of the quasilinear theory of the yl-effect 
with the results of the numerical simulations only positive 
values of the anisotropy a-parameter are reasonable. We shall 
see that for a > the rotation law in the supergranulation 
layer well approaches the observed state with dfl/dr < 
shown in Fig.^ The meridional flow at the surface is pole- 
ward and it is equatorward at the bottom of the convection 
zone. In particular, the resulting amplitude of the flow is 
an important outcome of the simulations. Doppler measure- 
ments by Komm et al. (1993 1 show a poleward flow on the 
surface, and the helioseismology indicates that this direction 
of the circulation persists to a depth of at least 12 Mm (Zhao 
& Kosovichev 2004 1. A characteristic amplitude for the sur- 
face velocity is 10 m/s. 

The global circulation computed with a = 2 is shown 
in the Figs.|5]and|6l The radial gradient of i? and the direc- 
tion of meridional flow in the subsurface region are close to 
the observed state. The computed rotation law reflects all the 
known basic features, i.e. 

- equatorial acceleration, 

- negative gradients of f? in the outer region, 

- small variation with depth for latitudes ^ 30°, 

- shallow tachocline 

(Fig. |5jl. The clear subrotation of the surface region 
idQ/dx < for x > 0.9) is mainly related to the ra- 
dial fluxes of angular momentum. Figure |3l shows that the 
functions Iq and Jq responsible for the radial fluxes are al- 
most completely quenched for fl* > 2, which corresponds to 
X < 0.9 (Fig. IS. Below this depth the results are thus unin- 
fluenced by the special choice of a. The subsurface increase 
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Fig. 5. Rotation law for a = 2. Left: angular velocity isolines. 
Right: depth profiles of the rotation rate. 
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Fig. 6. Streamlines of the meridional flow (left) and flow ve- 
locity for 45° latitude as function of depth (right) for a ~ 2. 
Positive velocity means poleward flow. One cell of anticlock- 
wise circulation occupies the entire convection zone. 

of rotation rate may be important for the solar dynamo (Bran- 
denburg 2005 1. 

The meridional flow at the surface is poleward. Its ampli- 
tude is slightly below 10 m/s. The circulation consists of a 
single cell with a return equatorward flow at the bottom. The 
flow reversal occurs rather deep so that the amplitude of the 
bottom velocity, ~10 m/s, is not small despite of the strong 
downward density increase. 

The poles are warmer than the equator by 2.3 K in the 
model given in Figs.|5land|6lfor a = 2. Many computations 
with various a-values were made. The models with a ~ 2 
reproduce the observations most closely. 

Angular momentum transport by rotating turbulence is 
known to originate from anisotropies existing in the turbu- 
lence and/or from any inhomogeneity of the turbulent fluid. 
For fast rotation the yl-effect produced by an inhomogeneity 
is much stronger compared to the effect of pure anisotropy. 
This is probably because the rotation can easily modify the 
anisotropies, but not the stratification. The modification of 
the spectral tensor introduced in this paper influences the 
anisotropy with consequences for the A-effect and the re- 
sulting rotation law close to the surface. By this modification 
the anisotropy parameter a enters the theory. The parameter 
is fixed using the results of 3D numerical box simulations. 
Global simulations have made considerable progress recently 
(Robinson & Chan l200l1 Brun & Toomre 120021 Brun, Mi- 
esch & Toomre '2004't but seem to suffer still from the pro- 
blem of resolution for the smaller scales. 
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Fig. 7. Depth of penetration of the meridional flow into the 
radiative zone at latitude 45° in dependence on the viscosity 
prescribed for the convectively stable region. The computed 
dependence is very close to Dpcn oc y^Iwc. 



equator temperature difference resulting from anisotropic dif- 
fusion inside the convection zone (Riidiger et al. '2005 a") is 
amplified. Our computations, however, show that the influ- 
ence of the region of stable stratification on the global flow in 
the convection zone reduces with decreasing i^coio and Xcoic^ ■ 
For the smallest i^coio given in Fig.Qthe differential rotation 
and the temperature distribution inside the convection zone 
are almost the same as in models without penetration. 

Our computations do not support any significant pene- 
tration of the meridional circulation beneath the convection 
zone. The calculations confirm our former result (Kitchatinov 
& Riidiger 1 995 ) that the penetration is very small unless the 
tachocline below the convection zone is strongly turbulent. 
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4.2. Penetration 

With the Model 2 the penetration depth of the meridional 
circulation into the radiative zone was computed. The pen- 
etrating flow shows multiple reversals of the meridional velo- 
city with increasing depth. The velocity amplitude decreases 
many times after each reversal. The penetration depth Dp^n 
of Fig.0is defined as the distance from the base of the con- 
vection zone to the location of the first reversal of the sign of 

Ug. 

The slope in the plot of Fig.0is very close to 0.5. This 
result complies with the finding of Oilman & Miesch ('2004!) 
that the penetration under solar conditions belongs to the Ek- 
man regime. The penetration of this type results from vis- 
cous drag imposed by the meridional flow at the base of the 
convection zone on the fluid below. The standard estimate 
for the Ekman depth is Dpcn ^ y^^torc/STi. This is further 
supported by our finding that the variation of heat conduc- 
tivity for constant ;/corc does not change Z3pon- The compu- 
tations, however, do not reproduce the expected dependence 
Dpen oc O^*^^. The slope is not constant and slightly larger 
than -0.5. The rotation rate dependence is better (though still 
not precisely) reproduced by the relation, 

-1/4 



(12) 



expected for nonuniform rotation. 

For large i^corc the penetration has a considerable effect 
on the global circulation and the thermodynamics inside the 
convection zone. The meridional flow at the base of the con- 
vection zone is reduced and the differential rotation and the 
differential temperature are increased with allowance for pen- 



etration. This is in agreement with Rempel ( I2005t who found 
that the balance of the meridional forces in the tachocline 
region demands a latitudinal entropy gradient which then 
spreads into the convection zone by radial diffusion. The 
poles become even warmer with this effect and the pole- 
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Appendix A: A-effect by rotating 
inhomogeneous turbulence 

Here we systematically derive the yl-effect effect for anelastic 
turbulent flow of a rotating stratified fluid. The usual strategy 
to deal with rotating turbulence is to prescribe the turbulence 
of the nonrotating fluid and then to derive the influence of 
rotation on the given 'original' turbulence. Here a nonrotating 
anelastic flow, div(ptt) = 0, is considered. 

The density stratification is basic for the existence of the 
vl-effect. Fourier modes of the momentum density, 

rh{z,uj) = J exp {—ix ■ z + iLut) 

■ m{x,t) da;dt/(27r)'', (Al) 
are used in order to construct the spectral tensor, 

{z,z',LU,u;') ^ {m, {z,Lu)rhj {z',lu')), (A2) 

of the fluctuating momentum density, m = pu. The turbu- 
lence is assumed statistically steady. The two new wave vec- 
tors, 



k = {z- z') /2, 
are introduced. Then 



K, = Z + Z 



Ivh, {z, z', Lj, iv') ^5{lj' + lj) Mlf (k, fc, UJ) 



(A3) 



(A4) 



The dependence on k describes the properties of the correla- 
tions on the small scale while the dependence on k accounts 
for the turbulence inhomogeneity on large scale (Roberts & 
Soward [T975 l). We further assume that the (nonhelical) turbu- 
lence does not possess any preferred directions other than its 
intensity gradient. The overall structure of the tensor, 



(0) 



1 



aSij + b—r- + c—TT (Kikj — Kjki) + d 



[k- k) 



(A5) 



differs from a former derivation (Kitchatinov I1987t by the 
last extra term. The symmetry condition M^^^ {k, fc, u) — 
Mj°^ {k, — fc, Lo) , is fulfilled. The divergence-free conditions 

ih + «;./2) M^°^ = 0, 
(fc,-«;,/2)M(°) = 0, 



(A6) 



lead to two vector equations with two nontrivial components 
each, hence 

(fc • k)^ 
k2 (fc . Kf 



0, 



a-2c + d 



^2 

6 + 2c + e— = 0, 
c-2d-2e = 0. 



(A7) 



The system is of rank 3 so that three unknowns can be ex- 
pressed in terms of the other two (say) c and e. It is conve- 
nient to express c and e in terms of the spectra E and Ei, 
i.e. 

E (k, fc, Lu) ~ El (k, k, uj) - ^) 



327rfc2 (1 



El (k, fc, u) 
^ ~ 167rfc2 
This leads to 



(A8) 



+ 



^- ^1 [2W - -21?- 




(A9) 



It is further assumed that E and Ei depend on the wave num- 
ber fc only rather than on the wave vector fc. The function E 
is the spectrum of turbulence intensity, 

00 00 

{m^ [x)) ^ {u^) ^ j J E{x,k,u}) dkdcu, (AlO) 


with 

E{x,k,Lu) = J exp {ix ■ k) E {k,, k,Lu) dK. (All) 

The spectrum Ei does not contribute to the turbulence inten- 
sity but it controls the anisotropy, (u^) / (ug). The dimension- 
less anisotropy parameter, a, can be written as 

J J Ei{x,k,uj)dkduj, (A12) 



(m^ (a;)) 



and it can depend on position. For the special case of 

- Ei/E = const = a/4, 

- i? cx (5 (fc — ^cm-) (mixing-length approximation), 

- V{nn?) ~ -~2g{m?) /L (g is radial unit vector and L is a 
constant length), 

the spectral tensor iA9\ can be integrated over wave space to 
give the one-point correlation tensor, i.e. 

(u^) /I 
= TTT2 — 772 o'^y + :f?f ^^n " mi. 



(A13) 



1 



{5ij - "igigj) . 



1+11,,JL^ 15 
For the rms intensities in the radial and the horizontal direc- 
tions it yields 



iul) = 



.1 A , 4orr2a 



lu^) 1 ^ ''^2 p2 



(A14) 
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APPENDIX A: A-EFFECT BY ROTATING INHOMOGENEOUS TURBULENCE 



The condition that the turbulence intensities (IA15> are 
positive-definite restricts the value of the anisotropy parame- 
ter a, 

15 ^^^^^ 
-27^<a<- + — . (A15) 

Isotropy ((u,^) = {ug)) results for a = 5/2. Although the 
turbulence is isotropic in this case, the rotation produces a 
vl-effect with negative Qr^- 

The ratio n/k in Eq. (IA9> is of order of icorr/L. For 
a weakly inhomogeneous turbulence, therfore, the terms of 
higher than second order in the scale ratio can be neglected 
so that 

Mio) _ E{K,k,uj) f _ A , A k^kJ 



Join*) 



i + n* 



n* 



arctan 



1 / , An* 



n*"^ - 1217*2-45 



n* 



1 + n*^ 

arctan f^* ) , (A18) 



and the expressions 

loin*) ^ 

X (-19- 



i + n* 



in*^ + 2A 
n* 



arctan n* , 



1/ 1 i\ f^it\jj 

El in, k,io) I ik ■ k) 



(A16) 
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+ 



ik ■ kY 



An*^ 

X (-15 



(A19) 



2n*^ 3^2*2 + 15 



arctanfJ* 



16^fc4 \^ fc2 r r.,, 

\2 ' 



:,2 



1 



ik-K Y 



i + n*^ n* 

are new. For slow rotation (il* ^ 1) only the radial flux of 
angular momentum exists, i.e. 



i^ij kikj/k 



^0-4 /o--^ Ji^h^o{n*^ 



(A20) 



results. This equation directly leads to the A-effect. After the 
standard derivations (KR93), one arrives at 



Note that Iq is negative for slow rotation. For fast rotation 
(n* ^ 1), Ji dominates all other functions. 



= z/T ^kgi {V^"'' in*) igiEjki + gjCiki) 

- in*) in,e,M + n,em) 



(A21) 



(A17) 



The azimuthal components of this tensor reproduce Eq. (|3|l 
with the normalized fluxes of angular momentum defined by 
Q and 0. The functions Jq and Ji Ql are the same as in 
KR93, 



The contribution of the anisotropy parameter a vanishes for 
fast rotation. 



